MSDS 6373 Mid Term Take Home Portion

URLs

R shiny : R shiny app
Markdown : R Markdown

Load Midterm2020 Dataset

msds_midterm_ds <- read.csv('data/midterm2020.csv')

Plot the data

ACF and Spectral Density using plotts.sample.wge Partial ACF

Condition 3 - Autocorrelation Function

acf(msds_midterm_ds$x[1:72],lag.max = 72)

acf(msds_midterm_ds$x[73:144],lag.max = 72)

Partial ACF shows lag 59 is non-significant that means 59th lag has some effect on current observation.

Stationarity

Condition 1: Mean does not depend on time. We don’t know what this dataset is about so cannot apply domain knowledge. But spectral density plot shows there are peaks at 0,0.08,0.18,0.25,0.33,0.43 and at 0.5 indicates seasonality present in the dataset which violates the condition of stationarity. It’s difficult to assess mean from such a small dataset but it is showing some upward trend for small duration for observations between 80-120 and the downward trend. So Mean may be depends on time.

Condition 2: Variance is finite and does not depend on time.

By looking at plot, it doesn’t look like it violates finite variance assumption. It looks like variance around mean is almost same across time so it may not be violating this assumption

Condition 3: Covariance depends only on how far apart in time not where in time.

ACF plots below almost looks same and there is very little evidence against constant covariance across time.

In general, ACF plot exhibits damped, oscillatory and sinusoidal behavior which is also sign of non-stationarity.

From visual evidence, due to seasonality and mean this time series does not appear to be stationary.

GLP form

Here are both models in GLP forms

Get the Psi weights and substitute in GLP equation

m1_phi   <- c(0.538,0.0606,0.1923)
m1_phi_1 <- mult.wge(fac1 =c(rep(0,11),1),fac2=m1_phi)
m1_psi   <- psi.weights.wge(phi=m1_phi_1$model.coef,lag.max = 4)

print(m1_psi)

## [1] 0.5380000 0.3500440 0.4132265 0.3469859

m2_phi   <- c(1.0507,-0.0756)
m2_theta <- c(0.5927,0.2751)
m2_psi   <- psi.weights.wge(phi=m2_phi,theta=m2_theta,lag.max = 4)

print(m2_psi)

## [1] 0.45800000 0.13052060 0.10251319 0.09784326

GLP equation

Model 1 :

\(X_{t}\) = \(a_{t} + 0.538 *a_{t-1} + 0.3500440 *a_{t-2} + 0.4132265*a_{t-3} + 0.3469859*a_{t-4} + ...+..\)

Model 2 :

\(X_{t}\) = \(a_{t} + 0.458*a_{t-1} + 0.1305206*a_{t-2} + 0.10251319*a_{t-3} + 0.09784326*a_{t-4} + ...+..\)

Model 2 Invertible ?

m2_phi   <- c(1.0507,-0.0756)
m2_theta <- c(0.5927,0.2751)

print(" AR portion ")

## [1] " AR portion "

factor.wge(phi=m2_phi)

## 
## Coefficients of Original polynomial:  
## 1.0507 -0.0756 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-0.9730B              1.0277               0.9730       0.0000
## 1-0.0777B              12.8704               0.0777       0.0000
##   
## 

print(" MA portion ")

## [1] " MA portion "

factor.wge(phi=m2_theta)

## 
## Coefficients of Original polynomial:  
## 0.5927 0.2751 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-0.8988B              1.1126               0.8988       0.0000
## 1+0.3061B             -3.2671               0.3061       0.5000
##   
## 

print("GLP equation upto first 4 terms  ")

## [1] "GLP equation upto first 4 terms  "

m2_psi   <- psi.weights.wge(phi=m2_phi,theta=m2_theta,lag.max = 4)
factor.wge(m2_psi)

## 
## Coefficients of Original polynomial:  
## 0.4580 0.1305 0.1025 0.0978 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-0.8812B              1.1349               0.8812       0.0000
## 1-0.0194B+0.2509B^2    0.0386+-1.9960i      0.5009       0.2469
## 1+0.4425B             -2.2598               0.4425       0.5000
##   
## 

AR and MA and ARMA Model is invertible as all roots are outside unit circle.

ACF and spectral Densities

Model 1

m1_phi  <- c(0.538,0.0606,0.1923)
model1r <- gen.aruma.wge(200,phi=m1_phi,d=0,s=12,sn=10)

plotts.sample.wge(model1r)

## $autplt
##  [1] 1.00000000 0.34685138 0.39281159 0.55900649 0.15809757 0.24293057
##  [7] 0.30029996 0.21899499 0.11746146 0.48440973 0.28784936 0.22892227
## [13] 0.83148323 0.23763288 0.31898189 0.47296387 0.08714661 0.16459815
## [19] 0.20100990 0.11864446 0.01580663 0.35413504 0.17430689 0.12361858
## [25] 0.68481365 0.15114099
## 
## $freq
##   [1] 0.005 0.010 0.015 0.020 0.025 0.030 0.035 0.040 0.045 0.050 0.055 0.060
##  [13] 0.065 0.070 0.075 0.080 0.085 0.090 0.095 0.100 0.105 0.110 0.115 0.120
##  [25] 0.125 0.130 0.135 0.140 0.145 0.150 0.155 0.160 0.165 0.170 0.175 0.180
##  [37] 0.185 0.190 0.195 0.200 0.205 0.210 0.215 0.220 0.225 0.230 0.235 0.240
##  [49] 0.245 0.250 0.255 0.260 0.265 0.270 0.275 0.280 0.285 0.290 0.295 0.300
##  [61] 0.305 0.310 0.315 0.320 0.325 0.330 0.335 0.340 0.345 0.350 0.355 0.360
##  [73] 0.365 0.370 0.375 0.380 0.385 0.390 0.395 0.400 0.405 0.410 0.415 0.420
##  [85] 0.425 0.430 0.435 0.440 0.445 0.450 0.455 0.460 0.465 0.470 0.475 0.480
##  [97] 0.485 0.490 0.495 0.500
## 
## $db
##   [1]  14.90482180   7.45126106  -2.17041089   0.58329601  -5.22051008
##   [6]   0.41124686  -4.79609505  -4.42875981  -9.01193150 -19.74796944
##  [11]  -4.86083851 -21.09774593 -11.34954015  -7.11760641  -8.15110351
##  [16]  10.47533036   5.27707819  -1.73023035  -7.01809263  -4.32122025
##  [21] -11.11492686  -8.46981404  -7.94306592  -9.19321070  -5.99941765
##  [26]  -9.61665589  -9.30683434  -8.19274013 -13.55280620  -9.36621935
##  [31]  -7.63188835  -9.65723532   4.09927882  -4.87140332 -17.23535514
##  [36]  -8.82213961 -18.73067921 -13.60772553 -15.87086658  -7.46272867
##  [41] -17.33235276 -12.77541423 -18.01382578 -10.13905891 -13.62032090
##  [46] -10.51538041 -15.00404268  -5.44990224   0.04660156   0.20393620
##  [51]  -1.67599919  -6.99300423 -15.49481823 -13.84948282 -10.22241165
##  [56] -18.46881119 -16.88096213 -20.89953000 -17.45210139 -16.93542192
##  [61] -13.93861343 -12.46526458  -5.59878442  -4.80820451  -2.87770207
##  [66]   5.63503108  11.57768742   0.17212329  -2.86679836  -8.88891394
##  [71] -11.67234863 -17.00978308 -10.65961079 -15.41762633 -24.56508809
##  [76] -21.16810572 -11.82379715 -20.74545137 -14.58017069 -22.70744258
##  [81] -13.13198890  -2.61235900   9.11565276   3.60161214  -1.85578434
##  [86] -10.06620493 -10.03428152 -12.16141595 -14.79079236 -18.14174711
##  [91] -18.49268285 -12.90474064 -25.80023907 -18.55365722 -15.64994986
##  [96] -13.40377797 -21.92733659 -22.45707710 -12.22062551   2.04382858
## 
## $dbz
##   [1]   8.60791113   8.25105199   7.65421899   6.81622646   5.74111470
##   [6]   4.44966729   3.00606345   1.56795012   0.42966334  -0.08307103
##  [11]   0.08585499   0.65984895   1.31992929   1.87542921   2.24196592
##  [16]   2.38684833   2.29711813   1.96591225   1.38797896   0.55993217
##  [21]  -0.51524668  -1.81639718  -3.27781641  -4.74202867  -5.92296614
##  [26]  -6.50892087  -6.43844952  -5.95885771  -5.36403849  -4.83329813
##  [31]  -4.44759827  -4.23989503  -4.22521600  -4.41299402  -4.81039097
##  [36]  -5.42108658  -6.23986197  -7.23917883  -8.33929185  -9.35722841
##  [41]  -9.98287127  -9.93210175  -9.23257092  -8.19248150  -7.10956119
##  [46]  -6.15095312  -5.38842063  -4.84755439  -4.53530037  -4.45266318
##  [51]  -4.60058943  -4.98219106  -5.59951443  -6.43470700  -7.38748577
##  [56]  -8.13278705  -8.05916787  -6.83441588  -4.92288663  -2.92808402
##  [61]  -1.14184147   0.34952919   1.53478667   2.42327431   3.02705834
##  [66]   3.35561702   3.41439840   3.20444856   2.72233479   1.96054926
##  [71]   0.90991072  -0.43081884  -2.02198485  -3.68750670  -4.94575894
##  [76]  -5.13497615  -4.23058384  -2.88061785  -1.57992455  -0.51328283
##  [81]   0.27591630   0.78623312   1.02378632   0.99243077   0.69149067
##  [86]   0.11555224  -0.74513461  -1.90335241  -3.37040561  -5.14095533
##  [91]  -7.14599177  -9.13398860 -10.53841343 -10.82234591 -10.21122969
##  [96]  -9.32826255  -8.54029466  -7.96506080  -7.62273987  -7.50964429

PACF

pacf(model1r,lag.max = length(model1r))

Factor table

m1_phi   <- c(0.538,0.0606,0.1923)
m1_phi_1 <- mult.wge(fac1=c(rep(0,11),1),fac2=m1_phi)
factor.wge(phi=m1_phi_1$model.coef)

## 
## Coefficients of Original polynomial:  
## 0.5380 0.0606 0.1923 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 1.0000 -0.5380 -0.0606 -0.1923 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-1.7321B+1.0000B^2    0.8660+-0.5000i      1.0000       0.0833
## 1-1.0000B+1.0000B^2    0.5000+-0.8660i      1.0000       0.1667
## 1+1.0000B             -1.0000               1.0000       0.5000
## 1+1.7321B+1.0000B^2   -0.8660+-0.5000i      1.0000       0.4167
## 1+0.0000B+1.0000B^2    0.0000+-1.0000i      1.0000       0.2500
## 1+1.0000B+1.0000B^2   -0.5000+-0.8660i      1.0000       0.3333
## 1-1.0000B              1.0000               1.0000       0.0000
## 1-0.8650B              1.1560               0.8650       0.0000
## 1+0.3270B+0.2223B^2   -0.7356+-1.9893i      0.4715       0.3064
##   
## 

Model 2

m2_phi   <- c(1.0507,-0.0756)
m2_theta <- c(0.5927,0.2751)
model2r <- gen.arma.wge(200,phi=m2_phi,theta=m2_theta,sn=10)

plotts.sample.wge(model2r)

## $autplt
##  [1] 1.00000000 0.57071300 0.24418912 0.18115183 0.20431659 0.27880340
##  [7] 0.25649078 0.16012999 0.13333694 0.18063709 0.17694172 0.16384519
## [13] 0.10418848 0.13285187 0.20357607 0.25124086 0.14019937 0.07192981
## [19] 0.05088910 0.12472428 0.22859891 0.19893306 0.15616702 0.13835504
## [25] 0.13217860 0.14664477
## 
## $freq
##   [1] 0.005 0.010 0.015 0.020 0.025 0.030 0.035 0.040 0.045 0.050 0.055 0.060
##  [13] 0.065 0.070 0.075 0.080 0.085 0.090 0.095 0.100 0.105 0.110 0.115 0.120
##  [25] 0.125 0.130 0.135 0.140 0.145 0.150 0.155 0.160 0.165 0.170 0.175 0.180
##  [37] 0.185 0.190 0.195 0.200 0.205 0.210 0.215 0.220 0.225 0.230 0.235 0.240
##  [49] 0.245 0.250 0.255 0.260 0.265 0.270 0.275 0.280 0.285 0.290 0.295 0.300
##  [61] 0.305 0.310 0.315 0.320 0.325 0.330 0.335 0.340 0.345 0.350 0.355 0.360
##  [73] 0.365 0.370 0.375 0.380 0.385 0.390 0.395 0.400 0.405 0.410 0.415 0.420
##  [85] 0.425 0.430 0.435 0.440 0.445 0.450 0.455 0.460 0.465 0.470 0.475 0.480
##  [97] 0.485 0.490 0.495 0.500
## 
## $db
##   [1]  13.67907098   0.43301725  -0.11356620  -1.46965499   5.87567782
##   [6]  -1.56437549  -3.91904709   6.25030887   1.61970472   3.05255137
##  [11]   3.39962384  -1.85702134  -0.63012228   3.00070715   3.65178834
##  [16]  -2.46929209 -14.77345451   1.04561203   3.82122721   0.72556170
##  [21]  -2.03903417  -0.58934802  -8.57977739 -12.94503406   0.63529090
##  [26]   4.07748026   3.99307152  -5.39391513   2.53473650   5.67632220
##  [31]   0.09952773 -17.91348374  -1.78627115  -4.81818426  -4.59554995
##  [36]  -4.92360349   3.87776532  -0.40759807   5.44875856   3.03326794
##  [41]  -0.79868945   3.35471670   2.27611174  -7.55581932 -11.65238617
##  [46]  -2.71614862  -1.73950822  -6.39927222  -0.65262607  -2.65067520
##  [51]  -3.41375183  -3.39460461  -3.09123794 -11.43675871  -2.18186273
##  [56]  -0.66461536  -6.38554621  -4.28939182 -19.75768222  -5.10972935
##  [61]  -5.65719090  -9.20341861   1.05733350 -13.82994128 -10.09778019
##  [66]  -8.84528020  -1.24118030  -8.88013124   1.42139197  -4.40052964
##  [71]  -4.13438714  -8.00018054  -2.79890588  -4.60515668 -17.60666632
##  [76] -13.28016075  -8.14860180  -2.10579864  -9.26630849  -2.20926151
##  [81] -11.10107659 -21.36427629  -6.98873329  -8.34681710  -3.65458391
##  [86] -15.06974361 -26.71919639  -4.91868755  -9.58590217 -13.17073430
##  [91]  -8.22741069   0.13134574  -6.06707604 -12.69980839  -7.75146344
##  [96]  -8.82988446 -10.64235290 -16.95356019 -14.14799613 -16.64327773
## 
## $dbz
##   [1]   7.40791448   7.15409095   6.74421789   6.20007229   5.55543948
##   [6]   4.85703391   4.16132925   3.52463243   2.98790762   2.56468199
##  [11]   2.24128883   1.98909915   1.77895573   1.58938649   1.40760828
##  [16]   1.22705970   1.04543191   0.86507487   0.69514514   0.55298268
##  [21]   0.46164012   0.44202290   0.50209229   0.62959140   0.79377786
##  [26]   0.95517339   1.07728120   1.13530004   1.12081961   1.04362144
##  [31]   0.93109452   0.82389664   0.76565023   0.78731393   0.89314155
##  [36]   1.05768530   1.23588151   1.37832528   1.44285425   1.39970070
##  [41]   1.23227818   0.93649541   0.52037885   0.00448354  -0.57759703
##  [46]  -1.18050596  -1.75282317  -2.24845166  -2.64016204  -2.92705226
##  [51]  -3.13019218  -3.27987659  -3.40324385  -3.51786525  -3.63129369
##  [56]  -3.74384283  -3.85193052  -3.95038934  -4.03332711  -4.09412160
##  [61]  -4.12568069  -4.12192384  -4.08054246  -4.00593131  -3.91060108
##  [66]  -3.81399811  -3.73906557  -3.70790875  -3.73789283  -3.83879453
##  [71]  -4.01103372  -4.24494323  -4.52142822  -4.81482237  -5.09857654
##  [76]  -5.35296123  -5.57169647  -5.76359601  -5.94771892  -6.14433271
##  [81]  -6.36561653  -6.60893695  -6.85420980  -7.06692457  -7.20846251
##  [86]  -7.25224139  -7.19799129  -7.07447850  -6.92913809  -6.81300890
##  [91]  -6.76977471  -6.83126234  -7.01699009  -7.33452169  -7.77819239
##  [96]  -8.32456723  -8.92409169  -9.49212948  -9.91177183 -10.06828950

PACF

pacf(model2r,lag.max=length(model2r))

Factor table

Factor table using Phis and thetas

m2_phi   <- c(1.0507,-0.0756)
m2_theta <- c(0.5927,0.2751)
factor.wge(phi=m2_phi)

## 
## Coefficients of Original polynomial:  
## 1.0507 -0.0756 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-0.9730B              1.0277               0.9730       0.0000
## 1-0.0777B              12.8704               0.0777       0.0000
##   
## 

factor.wge(phi=m2_theta)

## 
## Coefficients of Original polynomial:  
## 0.5927 0.2751 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-0.8988B              1.1126               0.8988       0.0000
## 1+0.3061B             -3.2671               0.3061       0.5000
##   
## 

Factor table using GLP equation upto first 4 terms

m2_phi   <- c(1.0507,-0.0756)
m2_theta <- c(0.5927,0.2751)
model2r <- psi.weights.wge(phi=m2_phi,theta=m2_theta,lag.max = 4)
factor.wge(phi=model2r)

## 
## Coefficients of Original polynomial:  
## 0.4580 0.1305 0.1025 0.0978 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1-0.8812B              1.1349               0.8812       0.0000
## 1-0.0194B+0.2509B^2    0.0386+-1.9960i      0.5009       0.2469
## 1+0.4425B             -2.2598               0.4425       0.5000
##   
## 

Average Squared Error

Model 1

This is AR(3) Model with s=12

m1_phi  <- c(0.538,0.0606,0.1923)
forecast_1 <- fore.aruma.wge(msds_midterm_ds$x,phi=m1_phi,s=12,n.ahead=12,lastn=TRUE,limits=F)

ASE1 <- mean((msds_midterm_ds$x[(144-12+1):144]-forecast_1$f)^2)
ASE1

## [1] 1301829

Model 2

This is ARMA(2,2)

m2_phi   <- c(1.0507,-0.0756)
m2_theta <- c(0.5927,0.2751)
forecast_2 <- fore.arma.wge(msds_midterm_ds$x,phi=m2_phi,theta=m2_theta,lastn=T,n.ahead=12)

ASE2 <- mean((msds_midterm_ds$x[(144-12+1):144]-forecast_2$f)^2)
ASE2

## [1] 836289.1

Rolling Window ASE

Model 1

phis = m1_phi
thetas = 0
s = 12
d = 0

trainingSize = 62
horizon = 12
ASEHolder1= numeric()

for( i in 1:(144-(trainingSize + horizon) + 1))
{
  
  forecasts = fore.aruma.wge(msds_midterm_ds$x[i:(i+(trainingSize-1))],phi = phis, 
                             theta = thetas, s = s, d = d,n.ahead = horizon,lastn=TRUE,
                             plot=FALSE )
  
  ASE = mean((msds_midterm_ds$x[(trainingSize+i):(trainingSize+ i + (horizon) - 1)] - forecasts$f)^2)
  ASEHolder1[i] = ASE

}

ASEHolder1

##  [1]  495169.60  484480.66  241370.06 1485990.20 1716994.82 1334663.43
##  [7] 1198688.37 1049539.00  901377.10  772935.40  743298.54  677725.94
## [13]  694641.45  713824.23  646789.81  219756.19  326635.46   93374.82
## [19]   50189.98   49408.56   48310.13   58180.43   61677.13   94348.31
## [25]  168108.13  176103.23  180276.05  126503.71   60543.93   70832.53
## [31]   64621.52   75438.72   71559.01   59317.32   52267.65   51445.73
## [37]   66938.85   71544.03  108617.57  197048.92  294072.83  313171.59
## [43]  184797.42  202770.67  187478.28  286481.65  246845.62  374324.07
## [49]  446966.18  497929.58  455726.59  403340.91  334810.17  433707.39
## [55]  433416.43  490755.25  547618.33  462331.95  433774.86  338255.36
## [61]  240323.14  228149.32  245338.50  342847.34  350809.48  389398.85
## [67]  407727.61  415883.15  504509.62  608431.57  719190.27

hist(ASEHolder1)

WindowedASE = mean(ASEHolder1)
WindowedASE

## [1] 388474.9

Model 2

phis = m2_phi
thetas = m2_theta
s = 0
d = 0

trainingSize = 62
horizon = 12
ASEHolder2 = numeric()

for( i in 1:(144-(trainingSize + horizon) + 1))
{
  
  forecasts = fore.aruma.wge(msds_midterm_ds$x[i:(i+(trainingSize-1))],phi = phis, 
                             theta = thetas, s = s, d = d,n.ahead = horizon,lastn=TRUE,
                             plot=FALSE )
  
  ASE = mean((msds_midterm_ds$x[(trainingSize+i):(trainingSize+ i + (horizon) - 1)] - forecasts$f)^2)
  ASEHolder2[i] = ASE

}

ASEHolder2

##  [1] 1000930.9 1017219.8 1044817.8 1422122.5 1297751.8 1311985.0 1766911.3
##  [8] 1344516.1 1251862.0 1125764.1 1055355.0  938136.1  839616.7  777339.6
## [15]  687821.0  592184.1  571635.8  550375.0  782855.6  521787.5  526857.2
## [22]  478213.7  468372.8  430850.1  418300.4  428264.1  451063.4  490196.7
## [29]  503264.1  408096.5  576232.8  402154.2  420859.4  433757.5  457060.6
## [36]  456170.8  494490.0  525277.3  550726.2  581472.1  628226.9  580250.9
## [43]  861769.0  643304.0  662008.3  695455.2  718437.1  781183.4  857882.4
## [50]  902396.9  935129.0 1029049.6 1041887.6 1023285.4 1349657.9 1106748.1
## [57] 1140600.8 1088221.9 1067370.1  914507.2  736515.2  625772.6  545726.1
## [64]  474438.5  465390.2  571170.0 1038769.7  789730.9  788529.9  782332.0
## [71]  835488.9

hist(ASEHolder2)

WindowedASE = mean(ASEHolder2)
WindowedASE

## [1] 775829.6

Window ASE comparison

library(plotly)

## Loading required package: ggplot2

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

x <- c(1:71)
data <- data.frame(x, ASEHolder1)
fig <- plot_ly(data, x = ~x, y = ~ASEHolder1, type = 'scatter', mode = 'lines', name = 'ARIMA(3,0,0) s =12')
fig <- fig %>% add_trace(y = ~ASEHolder2, name = 'ARMA(2,2)', mode = 'lines+markers')
fig <- fig %>% layout(title = "Window ASE comparison",
                      xaxis = list(title = "Iteration (Window)"),
                      yaxis = list (title = "Average Squared Error"))

fig